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ABSTRACT 

The energy and momentum deposited by the radiation from accretion flows onto the 
supermassive black holes (BHs) that reside at the centres of virtually all galaxies 
can halt or even reverse gas inflow, providing a natural mechanism for supermassive 
BHs to regulate their growth and to couple their properties to those of their host 
galaxies. However, it remains unclear whether this self-regulation occurs on the scale 
at which the BH is gravitationally dominant, on that of the stellar bulge, the galaxy, 
or that of the entire dark matter halo. To answer this question, we use self-consistent 
simulations of the co-evolution of the BH and galaxy populations that reproduce the 
observed correlations between the masses of the BHs and the properties of their host 
galaxies. We first confirm unambiguously that the BHs regulate their growth: the 
amount of energy that the BHs inject into their surroundings remains unchanged 
when the fraction of the accreted rest mass energy that is injected, is varied by four 
orders of magnitude. The BHs simply adjust their masses so as to inject the same 
amount of energy. We then use simulations with artificially reduced star formation 
rates to demonstrate explicitly that BH mass is not set by the stellar mass. Instead, 
we find that it is determined by the mass of the dark matter halo with a secondary 
dependence on the halo concentration, of the form that would be expected if the halo 
binding energy were the fundamental property that controls the mass of the BH. We 
predict that the black hole mass, mbh, scales with halo mass as TObh oc ?ri^ alo , with 
a w 1.55 ± 0.05 and that the scatter around the mean relation in part reflects the 
scatter in the halo concentration-mass relation. 

Key words: Cosmology: Theory - Galaxies: Active - Galaxies: Evolution - Galaxies: 
Formation - Hydrodynamics - Galaxies: Quasars: General 



1 INTRODUCTION 

Almost all massive galaxies are thought to contain a 
central supermassive black hole (BH) and the proper- 
ties of these BHs are tightly correlate d with those of 
the ga laxies in which the y resid e (e.g. Magorrian et al.l 
19981: iFerrarese fc Merrittl l200d: ICebhardt et al l l200d; 



Tremaine et all 120021 : lHaring fc RbJ 12004 iHopkins et al.l 
2007bl : iHd 12008 



It is known that most of the mass of 
the BHs is a ssembled via luminous accretion of matter 
l|Soltanl Il98l). The energy emitted by this process pro- 
vides a natural mechanism by which BHs can couple 
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e.g. ISilk & Reed Il998l: iHachnclt et al. Il998l: iFabian 


Adams et al 


l200ll: iKind 2003: Wvithe & Loeb 


iMurrav et al. 


20051: Merloni & Heinz! 2008), semi- 



analytic (e.g. iKauffmann fc Haehneltll200d ; ICattaned l200ll : 
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iGranato et al] 12004 I Bower et all l200d) and hydrody- 
namical (e.g. ISpringel et all l200d: iDi Matteo et all 120051: 



namical (e.g. [S pringcl et al. 2UUt>; Ui Ma tteo et al.l ZUUb; 
i Robertson et al.l l200d; ISiiacki et al.l 120071: I Hopkins et all 



2007al: IDi Matteo et all 120081 : lOkamoto et all 120081 : 
Booth fc Schavd 20091 ) studies have used this coupling 



between the energy emitted by luminous accretion and the 
gas local to the BH to investigate the origin of the observed 
correlation between BH and galaxy properties, and the 
buildup of the supermassive BH population. 

BHs are expected to regulate the rate at which they ac- 
crete gas down to the scale on which they are gravitationally 
dominant. For example, gas flowing in through an accretion 
disk can become so hot that its thermal emission becomes 
energetically important. Scattering of the photons emitted 
by the accreting matter by free electrons gives rise to the 
so-called Eddington limit. If the accretion rate exceeds this 
limit, which is inversely proportional to the assumed radia- 
tive efficiency of the accretion disk, then the radiative force 
exceeds the gravitational attraction of the BH and the in- 
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flow is quenched, at least within the region that is optically 
thin to the radiation. 

However, observations indicate that th e time-averaged 
accre tion rate is far below Eddington l|Kollmeier et al.l 
2006), suggesting the presence of processes acting on larger 
scales. Indeed, the existence of tight correlations between 
the mass of the BH and the properties of the stellar bulge 
indicates that sel f-regulation may happen on the scale of th e 
bulge (~ 1 kpc; I Adams et ail 1200 ll ; iHopkins et alj|2007a[ ), 
far exceeding the radius within which the BH is gravi- 
tationally dominant. However, since galaxy- wide processes 
such as galaxy mergers can trigger gas flows into the bulge 
l|Sanders et al.lll988l : [Mihos fc Hernauistl 19941 ). it is conceiv- 
able that BHs could regu late their growth on the scale of the 
entire galaxy (~ 10 k pc; lHaehnelt et~al1ll998l : iFabianll 19991 ; 
IWvithe fc Loebll2003l ) or even on that of the DM haloes host- 
ing t he galaxies (~ 10 2 kpc; I Silk fc Reesl 1 19981 ; iFerraresel 
|2002| ). Finally, it is possible, perhaps even likely, that self- 
regulation takes place simultaneously on multiple scales. For 
example, frequent, short, Eddington-limited outbursts may 
be able to regulate the inflow of gas on the scale of the bulge 
averaged over much longer time scales. 

In this paper we investigate, using self-consistent sim- 
ulations of the co-evolution of the BH and galaxy popula- 
tions, on what scale the self-regulation of BHs takes place. 
In Sec. [2] we describe the numerical techniques and simula- 
tion set employed in this study. In Sec. [3] we demonstrate 
that BH self-regulation takes place on the scale of the DM 
halo, and that the BH mass is determined by the binding 
energy of the DM halo rather than by the stellar mass of the 
host galaxy. Throughout we assume a flat ACDM cosmology 
with the cosmological parameters: {f2 m , S7b, J2a,<78,1s, h} = 
{0.238, 0.0418, 0.762, 0.74 0.951, 0.73|, as de termined from 
the WMAP 3-year data ijSpergel et al.ll2007h . 



2 NUMERICAL METHODS 

We have carried out a set of cosmological simulations us- 
ing Smoothed Particle Hydrodynamics (SPH) . We employ a 
significantly extended version of the parallel PMTr ee-SPH 
code GADGET ill (last described in ISpringe]||2005l ). a La- 
grangian code used to calculate gravitational and hydrody- 
namic forces on a particle by particle basis. The initial par- 
ticle positions and velocities are set at z = 127 using the 
Zeldovich approximation to linearly evolve positions from 
an initially glass-like state. 

In addition to hydrodynamic forces, we treat star for- 
mation, supernova feedback, radiative cooling, chemody- 
na mics and black hole accr e tion a n d feedback, as described 
in ISchave fc Dalla Vecchial (120081). iDalla Vecchia fc Schavel 
j200Sf).IWiersma et all (|2009al l. lWiersma et al.1 (|2009bl ) and 
iBooth fc Schavel (|20Qg| ) (hereafter BS09) respectively. For 
clarity we summarize here the essential features of the BH 
model, which i s itself a substantially modified version of that 
introduced bv lSprineel et al.1 J2005h . 

2.1 The black hole model 

Seed BHs of mass m scc( j = 10 _3 m g - where m g is the simu- 
lation gas particle mass - are placed into every DM halo that 
contains more than 100 DM particles and does not already 



contain a BH particle. Haloes are identified by regularly run- 
ning a friends-of-friends group finder on-the-fly during the 
simulation. After forming, BHs grow by two processes: accre- 
tion of ambient gas and mergers. Gas accretion occurs at the 
minimum of the Eddington rate, mEdd = 47rGmBHfnp/erCrTC 
and m accr = a47rG 2 m| H p/ (c 2 + v 2 ) 3 ^ 2 , where m p is the pro- 
ton mass, <tt is the Thomson cross-section, c is the speed of 
light, c s and p are the sound speed and density of the local 
medium, v is the velocity of the BH relative to the ambi- 
ent medium, and a is a dimensionless efficie ncy parameter. 
The p arameter a, which was set to 100 by ISpringel et al.l 
(2005), accounts for the fact that our simulations possess 
neither the necessary resolution nor the physics to accu- 
rately model accretion onto a BH on small scales. Note that 
for a = 1 this accretion rate reduces to the so called Bondi- 
Hoyle (|Bondi fc Hovldll944h rate. 

As long as we resolve the scales and physics relevant 
to Bondi-Hoyle accretion, we could set a — 1. If a simula- 
tion resolves the Jeans scales in the accreting gas then it 
will also resolve the scales relevant for Bondi-Hoyle accre- 
tion onto any BH larger than the simulation mass resolution 
(BS09). We therefore generally set a equal to unity. How- 
ever, this argument breaks down in the presence of a multi- 
phase interstellar medium, because our simulations do not 
resolve the properties of the cold, molecular phase, and as 
such the accretion rate may be orders of magnitude higher 
than the Bondi-Hoyle rate predicted by our simulations for 
star- forming gas. We therefore use a power-law scaling of the 
accretion efficiency such that a = (tih/tih)' 3 in star-forming 
gas, where = 0.1cm -3 is the critical density for the for- 
mation of a cold, star-forming gas phase. The parameter /3 
is a free parameter in our simulations. We set /3 = 2, but 
note that the results shown here are insensitive to changes 
in this parameter when /3 > 2 (see BS09), because in that 
case the growth of the BHs is limited by feedback. 

Energy feedback is implemented by allowing BHs to in- 
ject a fixed fraction of the rest mass energy of the gas they 
accrete into the surrounding medium. The energy deposition 
rate is given by 

lip ■ 2 e f e r .2 /,n 

b, = £f£ r m accr c = m B HC , (1) 

1 6 r 

where e r is the radiative efficiency of the BH, rh acC r is the 
rate at which the BH is accreting gas, and mBH is the rate 
of BH mass growth. 

We set e r to be 0.1, the mean value for ra- 
diatively efficient accret ion onto a Schwarzschild BH 
jShakura fc Sunvaevl l 1973b . We vary a but use et = 0.15 as 
our fiducial value. It was shown in BS09 that, for ei = 0.15, 
simulations identical to these reproduce the observed red- 
shift zero ?tibh — rn* and mBH — o relations, where a is the 
one-dimensional velocity dispersion of the stars and m* is 
the galaxy stellar mass. Energy is returned to the surround- 
ings of the BH 'thermally', that is, by increasing the tem- 
perature of A^oat of the BH's neighbouring SPH particles by 
at least AT m i n . A BH performs no heating until it has built 
up enough of an energy reservoir to heat by this amount. 
The use of an energy reservoir is necessary in these simula- 
tions as otherwise gas will be able to radiate away the energy 
every timestep. Imposing a minimum temperature increase 
ensures that the radiative cooling time is sufficiently long 
for the feedback to be effective. In our fiducial model we set 
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Table 1. Numerical parameters of the simulations. From left to 
right: Simulation identifier, comoving box size (Mpc/?t), number 
of both gas and DM particles, final redshift, gas particle mass 
(10 7 M Q //i), DM particle mass (10 7 Mq//i), maximum physical 
gravitational softening (kpc//i). Each simulation was run multiple 
times using different values of ef. 



Name 


L 




H 


m s 




€ max,phys 


L050N256 


50.0 


256 3 





8.7 


41.0 


2.0 


L 05 ON 128 


50.0 


128 3 





69.6 


328.0 


4.0 


L012N256 


12.5 


256 3 


2 


0.1 


0.6 


0.5 



ATheat = 1 and AT m i n = 10 K but the results are insensitive 
to the exact values of these parameters (see BS09). 



d 



x 
x 




2.2 The simulation set 

The simulations employed in the current work use cubic 
boxes of size 12.5 and 50 comoving Mpc//i and assume 
periodic boundary conditions. Each simulation contains ei- 
ther 128 3 or 256 3 particles of both gas and collisionless cold 
DM. Comoving gravitational softenings are set to 1/25 of 
the mean interparticle separation down to z = 2.91, below 
which we switch to a fixed proper scale. The 12.5 Mpc/h 
(50 Mpc/h) boxes are evolved as far as redshift two (zero). 
The numerical parameters of the simulations used in this 
study are summarized in Table [T] All results presented in 
this letter are derived from the 50.0 Mpc//i, 256 3 particle 
simulations, with the other box sizes and particle numbers 
employed to demonstrate numerical convergence. 



3 RESULTS AND DISCUSSION 

It is instructive to first consider under what conditions BHs 
can regulate their growth. To regulate its growth on a mass 
scale M sr , a BH of mass ttibh must be able to inject en- 
ergy (or momentum) at a rate that is sufficient to counter- 
act the force of gravity on the scale M SI , averaged over the 
dynamical time associated with this scale. The mass M SI 
could, for example, correspond to that of the BH, the stel- 
lar bulge, or the dark matter (DM) halo. If the BH cannot 
inject energy sufficiently rapidly, then gravity will win and 
its mass will increase. Provided that the maximum rate at 
which it can inject energy increases with ttibh (as is for ex- 
ample the case for Bondi-Hoyle and Eddington-limited ac- 
cretion with a constant radiative efficiency) and provided 
that this rate increases sufficiently rapidly to counteract the 
growth of M ar , the BH will ultimately reach the critical mass 
m-BH,crit(Af sr ) required to halt the inflow on the scale M sr . If, 
on the other hand, mBH S> m-BH,crit(M sr ), then the BH will 
quickly quench the accretion flow and its mass will conse- 
quently remain nearly unchanged. The BH will in that case 
return to the equilibrium value mBH,ciit(M SI ) on the time 
scale which characterises the growth of M ar . 

If the BH regulates its growth on the mass scale M SI 
and if moH <C M sr , then the critical rate of energy injection 
required for self-regulation is independent of the mass of 
the BH. It then follows from Eq. Q]that ttibh oc e^ 1 , which 
implies 
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Figure 1. Predicted redshift zero global BH mass density (black 
diamonds) and normalization of the mBH — a relation (black plus 
signs) as a function of the assumed efficiency of BH feedback, ef . 
Both quantities are normalized to their values in the simulation 
with ef = 0.15, which reproduces the observed relations between 
the mass of the BH and properties of the stellar bulge. Each 
point represents a different simulation. For 10 -4 < ef < 1 all 
data points track the dotted black line, which is a power-law 
with index minus one. This implies that in this regime BH mass 
is inversely proportional to rriBH; ar >d thus that the BHs inject 
energy into their surroundings at a rate that is independent of ef , 
as expected for self-regulated growth on scales that are sufficiently 
large for the gravity of the BH to be unimportant. The red data 
points show results from simulations with a mass resolution that 
is 8 times worse than the fiducial simulation. The blue data points 
correspond to simulations with 64 times better resolution than our 
fiducial resolution, but show results for redshift 2 rather than zero. 
The agreement between the black, red and blue points confirms 
numerical convergence and demonstrates that the BHs are already 
self-regulating at redshift 2. 

where m SOG d is the initial mass of the BH. Hence, if the self- 
gravity of BHs is negligible on the maximum scale on which 
they regulate their growth and if ttzbh ^> rn sco( i, then we 

expect TTlBH OC €f. 

The black diamonds plotted in Fig. [1] show the pre- 
dicted global mass density in BHs at redshift zero as a func- 
tion of ef, the efficiency with which BHs couple energy into 
the ISM, normalised to the density obtained for ef = 0.15. 
Similarly, the black plus signs indicate the normalisation of 
the m,BH — a relation divided by that for the ef = 0.15 run. 
The feedback efficiency, ef, is varied, in factors of 4, from 
ef = 9.2 x 10 -6 to ef = 9.6, which implies that the frac- 
tion of the accreted rest mass energy that is injected (i.e. 
e r ef) varies from 9.2 x 10~ 7 to 0.96. BH mass is clearly in- 
versely proportional to the assumed feedback efficiency for 
10 -4 < ef < 1. For ef > 1 the trend breaks down because 
the BH masses remain similar to the assumed seed mass, 
in accord with Eq. [2] If we had used a lower seed mass, 
then the trend would have extended to greater values of ef. 
The deviation from inverse proportionality that sets in be- 
low ef = 10 -4 is more interesting. Such low values yield BH 
masses that are more than 0.15/10 -4 ~ 10 3 times greater 
than observed, in which case they are no longer negligible 
compared to the masses of their host galaxies. In that case 
the critical rate of energy deposition will no longer be inde- 
pendent of mBH and we do not expect Eq. [2] to hold. 

We have thus confirmed that feedback enables BHs to 
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Figure 2. Median r«BH — "ihalo (left panel) and meH — m* (right panel) relations for all BHs more massive than 10m see ^. The black 
curves correspond to a simulation using our fiducial star formation law and the red, dashed curves show the result for a run in which the 
star formation efficiency was decreased by a factor of 100. In order to isolate the effect of stellar mass, we turned off supernova feedback in 
both runs. The BH scaling relations therefore differ somewhat from those predicted by our fiducial model, which does include supernova 
feedback. Baryons dominate the gravitational potential in the central regions of the galaxy when we use our fiducial star formation law, 
but DM dominates everywhere in the run with the reduced star formation efficiency. While the ?tibh — "i* relation is strongly affected 
by the change in the star formation efficiency, the relation between BH and halo mass remains invariant. This demonstrates that the BH 
mass is insensitive to the mass distribution on scales where the stellar mass dominates, and must instead be determined by the mass 
distribution on larger (^> lOkpc) scales. 



regulate their growth. Moreover, we demonstrated that this 
self-regulation takes places on scales over which the grav- 
itational influence of the BHs is negligible, provided that 
the fraction of the accreted rest mass energy that is coupled 
back into the interstellar medium is > 10 . 

To test whether it is the stellar or the dark matter distri- 
bution that determines the mass of BHs, we compare the BH 
masses in two simulations that are identical except for the 
assumed efficiency of star formation. One uses our fiducial 
star formation law, but in the other simulation we reduced 
its amplitude by a factor of 100, making the gas consump- 
tion time scale much longer than the age of the Universe. 
Because changing the amount of stars would imply chang- 
ing the rate of injection of supernova energy, which could 
affect the efficiency of BH feedback, we neglected feedback 
from star formation in both runs. In the simulation with 
'normal' star formation the central regions of the galaxies 
are dominated gravitationally by the baryonic component 
of the galaxy, whereas in the simulation with reduced star 
formation the DM dominates everywhere. Fig. [2] shows the 
m-BH — TTihaio and dibh — m* relations at redshift 0. While 
the two runs produce nearly identical BH masses for a fixed 
halo mass, the meH — m, relation is shifted to lower stellar 
masses by more than an order of magnitude in the model 
with reduced star formation. The insensitivity of the rela- 
tion between mBH and mhaio to the assumed star formation 
efficiency demonstrates that the BH mass is not set by the 
gravitational potential on the scale of the galaxy. We have 
verified that the same result holds at redshift two for the 
simulations with 64 times better mass resolution. Clearly, 
stellar mass does not significantly influence the relation be- 
tween the mass of the BH and that of its host halo. This 
implies that BH self-regulation occurs on the scale of DM 
haloes. 

If the rate by which the BHs inject energy is indepen- 
dent of the assumed feedback efficiency, then we expect the 
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Figure 3. The relation between BH mass and DM halo mass for 
all BHs that belong to central galaxies and have masses greater 
than 10m sco( j. The DM halo mass, rri200> is defined as the mass 
enclosed within a sphere, centred on the potential minimum of the 
DM halo, that has a mean internal density of 200 times the critical 
density of the Universe. The grey pixels show the results from our 
fiducial simulation (ef = 0.15), with the colour of each pixel set by 
the logarithm of the number of BHs in that pixel. The solid, red 
line shows the observational determination of the ttibh ~ "ihalo 
relation (Bandara et al. 2009) and has a slope of 1.55. The dotted, 
red lines show the la errors on the observations. The simulation 
agrees very well with the observed relation. The value of the slope 
and the scatter (which correlates with the concentration of the 
DM halo) suggest that the halo binding energy, rather than mass, 
determines the masses of BHs. 

same to be true for the factor by which BH feedback sup- 
presses star formation. This is confirmed by comparison of 
the global SFRs in runs with different values of ef (see Fig. 6 
of BS09). 

Fig. [3] compares the predicted log ln m B H - log 
relation with observation l|Bandara et al.|[2009l ^ . The agree- 
ment is striking. The slope and normalization of the ob- 
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served log lo (m,BH/M0) — log 10 (m ha io/10 13 M Q ) relation are 
1.55 ± 0.31 and 8.18 ± 0.11 respectively, whereas the sim- 
ulation predicts 1.55 ± 0.05 and 8.01 ± 0.04. Note that the 
simulation was only tuned to match the normalization of the 
relations between tobh and the galaxy stellar properties. 

If the energy injected by a BH is proportional to the 
halo gravitational bin ding energy, then, for isothermal mod- 
els (jSilk fc Rees!ll998h . rriBH oc mJ , . Here we extend these 
models to the more realistic universal halo density profile 
l|Navarro et al.l Il997l ). whose shape is specified by a con- 
centration parameter, c (we assumed c oc ii^/nj, where 
Wmax an d «v are the maximum halo circular velocity and 
the circular velocity at the virial radius respectively). It 
is known that conc entration decreases with increasing hal o 
mass, c oc m^ 1 (|Bullock et all l200ll ; Ipuffv et al.ll2008h , 
which then affects BH mass through the dependence of halo 
binding energy on concentration. If the total energy injected 
by a BH of a given mass is proportional to the energy re- 
quire d to unbind gas from a DM halo l|Lokas fc Mamonl 
l200ll ) out to some fraction of the virial radius, r j/r v then 



TTlBH oc 



(ln(l + c) -c/(l + c)) 

1 21n(l + c^) 



1 - 



(l + c, 



n 5/3 



(3) 



Inserting c oc m~ ' and computing the logarithmic deriva- 
tive with respect to m v in the mass range 10 10 Mq < m v < 
10 14 Mq, we find that the slope is a weak function of r c j/r v 



that varies from 1.50 at r cj 



The close match between theory, simulation and observa- 
tion suggests that the halo binding energy, rather than halo 
mass, determines the mass of the BH. 

The residuals from the dibh — w-haio relation 
(Alo g 10 771bh) are correlated with halo concentration 
(Spearman rank correlation coefficient p = 0.29, probability 
of significance P = 0.9998) as would be expected if 771bh 
is sensitive to the halo binding energy. The residuals are 
also correlated with galaxy stellar mass, though much less 
strongly (p = 0.09; P = 0.96). Taken together, these corre- 
lations tell us that, at a given halo mass, galaxies with BHs 
more massive than the average will also contain a larger 
than average amount of stars, and are hosted by more con- 
centrated haloes. This suggests that the galaxy stellar mass 
is also determined by the halo binding energy. Thus, outliers 
in the mBH — m-haio relation may still lie close to the mean 
wibh — m-j, relation. Furthermore, higher concentrations im- 
ply earlier formation times and spheroidal components do 
indeed typically host old stellar populations. 

In addition to the 'quasar mode' of feedback discussed 
in this work, it has recently become clear that a second 
'radio mode' may be required to qu ench cooling flows in 
galaxy groups and clusters (see e.g. ICattaneo et alj 120091 . 
for a review). Although we do not explicitly include a 'radio 
mode' in the current work, the AGN feedback prescription 
explored here is capable of suppressing cooling flows, at least 
on group scales, providing excellent matches to observed 
group density and temperature pr ofiles as well as galax y 
stellar masses and age distributions ^McCarthy et al . 2009). 
It is kn own that BH s obtain most of their mass in the 'quasar 
mode an 1982) so any discussion of what detemines the 
masses of BHs must focus primarily on this mode of accre- 



tion. Finally, the ability of a BH to quench cooling flows in 
the 'radio mode' is expect ed to be closely relate d the virial 
properties of the hot halo (jCattaneo et al.ll2009l ) and would 
therefore provide an additional link between BHs and DM 
haloes over and above what we discuss here and so serve 
to make any fundamental connection between BH mass and 
the properties of the DM halo even stronger. 

We conclude that our simulation results suggest that 
in order to effectively halt BH (and galaxy) growth, gas 
must not return to the galaxy on a short timescale. This 
requires that the BH injects enough energy to eject gas out 
to scales where the DM halo potential is dominant. The 
mass of the BH is therefore determined primarily by the 
mass of the DM halo with a secondary dependence on halo 
concentration, of the form that would be expected if the BH 
mass were controlled by the halo binding energy. The tight 
correlation between ttibh and m, is then a consequence of 
the more fundamental relations between halo binding energy 
and both 771bh and m». 
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